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The gap generation is studied in suspended clean graphene in the continuum model for quasipar- 
ticles with the Coulomb interaction. We solve the gap equation with the dynamical polarization 
function and show that, comparing to the case of the static polarization function, the critical cou- 
pling constant lowers to the value a c = 0.92, which is close to that obtained in lattice Monte Carlo 
simulations. It is argued that additional short-range four-fermion interactions should be included 
in the continuum model to account for the lattice simulation results. We obtain the critical line 
in the plane of electromagnetic and four-fermion coupling constants and find a second order phase 
transition separating zero gap and gapped phases with critical exponents close to those found in 
lattice calculations. 



I. INTRODUCTION 

Graphene, a one-atom-thick layer of graphite, is a remarkable system with many unusualproperties that was 
fabricated for the first time five years ago Theoretically it was shown long time ago Q that quasiparticle 
excitations in graphene have a linear dispersion at low energies and are described by the masslcss Dirac equation in 
2+1 dimensions. The observation of anomalous integer quantum Hall effect in graphene Q is in perfect agreement with 
the theoretical predictions Q and became a direct experimental proof of the existence of gapless Dirac quasiparticles 
in graphene. 

The unusual band structure of graphene has an important consequence for the electron-electron interaction in this 
material. In the continuum limit, graphene model on a honeycomb lattice, with both on-site and nearest-neighbor 
repulsions, maps onto a (2 + l)-dimensional field theory of Dirac fermions interacting through the Coulomb potential 
plus, in general, some residual short-range interactions represented by local four-fermion terms. The vanishing density 
of states at the Dirac points ensures that the Coulomb interaction between the electrons in graphene retains its long- 
range character due to vanishing of the static polarization function for q — > [|| . The large value of the "fine structure" 
coupling constant a = e 2 /hvp ~ 1 means that a strong attraction takes place between electrons and holes in graphene 
at the Dirac points. As is known, for graphene on a substrate with the effective coupling a/ft <C 1, k being a dielectric 
constant, the system is in a weak coupling regime and exhibits semimetallic properties due to the absence of a gap 
in the electronic spectrum. Much less is known about suspended graphene where the coupling constant is large. In 
fact, suspended graphene provides a condensed-matter an alog ue of strongly coupled quantum electrodynamics (QED) 
intensively studied in the 70-ties and 80-ties [1, 0, H, H, Hoj - The dynamics of the vacuum in QED leads to many 
peculiar effects not yet observed in nature. Some QED-like effects such as zitterbewegung (trembling motion) [111 ] . 
Klein tunneling [12] , Schwinger pair production [r| , supercritical atomic collapse fl4[_[l5| , have a chance to be tested 
in graphene (for experimental observation of the Klein tunneling in graphene, seejlqjv To observe these effects in 
graphene, it is important to use suspended and clean samples where charges from a substrate do not interfere with 
the dynamics of electrons. 

Recently, it is was shown in Ref. [13] that, for strong enough coupling a > a c , there is a tachyonic solution in the 
spectrum of the Bethe-Salpeter (BS) equation for the electron-hole bound state signaling the presence of excitonic 
instability of the zero-gap ground state of monolayer graphene in the supercritical regime. The critical coupling equals 
a c = 1/2 if the vacuum polarization is neglected and a c 1.62 in the random phase approximation with the static 
polarization [l8[. It was also shown there that physically the excitonic instability is connected with the well-known 
supercritical Coulomb center problem fl9l where a c = 1/2 in two spatial dimensions (20j . The situation is similar to 
that in the theory of superconductivity [21[ , where the four-fermion vertex instability has its origin in the Cooper pair 
problem. It was argued in [l7j that the formation of an excitonic condensate of electron-hole pairs should cure the 
excitonic instability and lead to opening of a quasiparticle gap in a free standing clean graphene resulting in dramatic 
changes in the transport properties. A similar situation occurs in QED in 3+1 dimensions where the gap generation 
takes place in the strong coupling regime H, 0, E3] (see, also, [1, Q). 

The problem of gap generation in graphene was considered before the actual fabrication of this material in Rcfs. 
[22l . I23I [24| where the random phase approximation with the static polarization function was used. Recently, lattice 
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Monte Carlo simulations found the value of the critical coupling a c = 1.08 for a semimetal- insulator transition [25[ 
and this motivated us to reconsider the problem of gap generation in graphene. It was already indicated in [23j that 
taking into account the frequency dependent polarization function should lower the critical coupling. We investigate 
this question in the present paper and confirm that the dynamical polarization is indeed quantitatively important 
: solving the gap equation with frequency dependent polarization function we find the critical coupling a c = 0.92 
instead of a c = 1.62 in the case of static polarization. We would like to note also that the presence of a gap would be 
valuable for electronics applications, in particular, for working graphene transistors. 

Another problem studied in the present paper is the order of phase transition connected with the gap generation 
in graphene. Due to the scale invariance of the model with the Coulomb interaction an infinite order phase transition 
was found in [23|, [24| ■ Such a phase transition belongs to the class of the so-called conformal phase transitions [2(| . 
According to the recent Monte Carlo (MC) simulations [2{| (for related MC simulations, see Ref . p7| ) . the semimetal- 
insulator phase transition in graphene is of the second order. One of the reasons for such a difference might be 
lattice finite size effects which can change the order of phase transition [28|, [2^, [3(| ■ On the other hand, according 
to [HI, HI, , the effective continuum theory for quasiparticles in graphene should contain besides the Coulomb 
interaction some additional contact four-fermion interaction terms that arise from the microscopic graphene lattice 
interactions. These terms contain a dimensionful parameter, therefore, they explicitly break the scale invariance of the 
continuum model. In such a case, one may expect a conventional second order phase transition. In order to take into 
account these four-fermion interaction terms, we consider in the present paper the simplest Gross-Neveu interaction 
term and show that the presence of this interaction term plays an important role: First, instead of a critical point 
we now have a critical line in the plane of electromagnetic and four-fermion coupling constants separating symmetric 
and symmetry broken phases. Second, the inclusion of this term indeed changes the order of phase transition from 
infinite to the second order along a part of the critical line < a < a c . Third, it lowers the value of the critical 
electromagnetic coupling comparing to the case of purely Coulomb interaction. At last, the critical indices stay closer 
to those obtained in lattice simulations (25[. 

The structure of the paper is the following. We begin with presentation in Sec. |TT]of the continuum model describing 
graphene quasiparticles interacting through the Coulomb potential. In Sec. IIHI we solve the gap equation with the 
frequency dependent one-loop polarization function and determine a critical coupling for the onset of a gap. To get 
insight into analytical solutions of the gap equation, we then turn back to the case of the static polarization and 
find asymptotical behavior of the gap function, calculate the excitonic condensate ("I"!') of particle-hole pairs, the 
correlation length, and critical exponents near the phase transition point a c . In Sec. [V]we include the Gross-Neveu 
four-fermion interaction, find explicitly the critical line in the plane of Coulomb and four-fermion interaction coupling 
constants, and determine the critical exponents for the phase transition along this line. In Conclusion we summarize 
the main results. 



II. THE MODEL 

For the description of the dynamics in graphene, we will use the same model as in Refs.[22l. [23j | in which while 
quasiparticles are confined to a two-dimensional plane, the electromagnetic (Coulomb) interaction between them is 
three-dimensional in nature. The low-energy quasiparticles excitations in graphene arc conveniently described in 
terms of a four-component Dirac spinor ty^ = (ipKAajipKBajipK'BajipK'Aa) which combines the Bloch states with 
spin indices a = 1, 2 on the two different sublattices (A, B) of the hexagonal graphene lattice and with momenta near 
the two-nonequivalent valley points (K,K') of the two-dimensional Brillouin zone. In what follows we treat the spin 
index as a "flavor" index with Nf components, a = 1, 2, . . . Nf, then Nf = 2 corresponds to graphene monolayer while 
Nf — 4 is related to the case of two decoupled graphene layers, interacting solely via the Coulomb interaction. 

The action describing graphene quasiparticles interacting through the Coulomb potential has the form 



dtd 2 r^ a (t, r) (ij°d t - WftV) r) 
- i J dtdt'd 2 rd 2 r'^ a {t 1 r) 7 °* a (t, r)U (t - t', |r - v'\)%(t', r')j°^ b (t', r'), (2.1) 

where vp is the Fermi velocity, \[r — vl/ty^ anc [ the 4x4 Dirac 7-matrices 7^ = r 3 (^) (a 3 , ia 2 , — ia 1 ) furnish a reducible 
representation of the Dirac algebra in 2 + 1 dimensions. The bare Coulomb potential Uo(t, |r|) is given by 

, , „ e 2 S(t) fd 2 ke^ r e 2 S(t) , . 

^W) = ^ /_- r = - j ^. (2.2) 
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However, the polarization effects considerably modify this bare Coulomb potential and the interaction will be 

U{t ^ r \>- K J 2 irJ 2tt |k|+n( W ,k) ' [Z - 6) 

where n is the dielectric constant due to a substrate and the polarization function H(uj, k) is proportional (within the 
factor 2tt/k) to the time component of the photon polarization function. Correspondingly, the Coulomb propagator 
has the form 

"<-•*" R+BW (2 ' 4) 

where the one-loop polarization function is [|| 

x ne 2 N f k 2 , . 

and in the instantaneous approximation it becomes 

H(. = 0,k) = ^|k|. (2.6) 

The continuum effective theory described by the action (|2.1[) possesses U{2Nf) symmetry. In the case of graphene, 
Nf = 2, the corresponding 16 generators are (see, for example, Ref . [23| ) : 

a a a a a a <j a 

—— ® I 4 , — ® 7 3 , ^r®7 5 , -5-877, (2-7) 

where I4 is the 4x4 Dirac unit matrix, and a a , with a = 0, 1, 2, 3, are four Pauli matrices connected with the spin 
degrees of freedom (er is the 2x2 unit matrix). However, as was pointed out in Ref. [3l| (see also Refs. (32l. l33j). 
this symmetry is not exact in the graphene tight-binding model on lattice. In fact, there are small on-site interaction 
terms which break the U(2Nf) - symmetry, their role will be considered in Sec. [V] 

III. GAP GENERATION AND THE CRITICAL COUPLING CONSTANT 

In this section we study spontaneous generation of a gap in the quasiparticle spectrum of graphene. The Schwinger- 
Dyson equation for the quasiparticle propagator has the form, 

/d 3 k 
j^D( P0 - k , p - k) 7 °S(fco, k)7°, (3.1) 

where the Coulomb propagator D(qq, q) is given by Eq. (|2.4p and in the random phase approximation the polarization 
is taken as in Eq. (|2.5[) . The vertex corrections are rather small § and we neglect them in what follows. 
The general form of the propagator of quasiparticles is 

5- 1 (p ,p) = ^ 1 Po7°-^P7-A, (3.2) 

where Z, A, A are functions of po,p and we included also a bare gap Ao. We assume that a dependence of these 
functions on the energy po is rather weak so that we can approximate these functions by their values at po = 0. In 
this approximation it is easy to see that Z = 1, then after the Wick rotation, fco = iw, we get a coupled system of 
equations for A(p), A(p): 

on 

r d 2 k . , . pk^(k) . . 

7^(^P- k ) ^ + k2 ^ (k)+AW ( 3 - 3 ) 
A(p) = A +'— I doj J ^^.p-k) — (3 . 4) 
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We write the integral over uj as 

00 oo 

f If dxf(x) \/x 2 + 1 

1 = / divDiuj, q)— 5 — t-^ = / -5—5 — , „ .„ T-g-, /(a;) = —. = , g = irN f a/4. (3.5) 

7 v ; w 2 + k 2 A 2 + A 2 y x 2 q 2 + k 2 A 2 + A 2 v 7 ^z 2 + 1 + y J ' v ; 

—00 —00 

The function f(x) changes slowly from 1/(1 + g) at x = (the instantaneous approximation for D(u>, q)) up to 1 at 
x = 00. The integral in Eq. (|3.5p can be evaluated exactly, 

T 1 \ J Vk 2 ^ 2 + A 2 (d 2 -l)(n~gc(d))+dg 2 c(g) 

1 = WWWTW J{d ^ d = — h\ — ' J{d ' g) = d^f-i ' (3 ' 6) 



. „ 2 cosh 1 (x) , . 2 cos l (x) , . , 

c(x) = V 2 , z>l, c(g) = L i, x<l, c(l) = 2. (3.7) 

vr — 1 VI — x 



where 



For A = and setting A = 1 on the right hand side of Eq. (|3.3j) we get the leading one- loop correction |34( , which 
comes from the range of momenta k p in the integral, 

= 1 + -T^T i 7T ~ 2 9+ (S 2 - !) c (5)) In - + finite terms, (3.8) 

where A is a momentum cutoff of order the inverse lattice spacing in graphene. The function A(p) renormalizes the 
Fermi velocity v* F {p) = vpA(p). The growth of v* F {p) in the infrared stops when a non-zero quasiparticle gap is taken 
into account (see, Eq. (|3.3[l ). In what follows we assume that the velocity renormalization is already performed (35| 
and put A = 1 in Eq. (|3.4p which then takes the form 

. . . . e 2 f d 2 k A(k) r/ , Ikl . ,,. 

Ap = A + - — 2- Y Jd= II, g , 3.9 

kJ (2tt) 2 | p -k| v / k 2 + A 2 (k) |p-k| 

where we set also A = in the variable d. Since the function J depends weakly on the angle between the vectors p 
and k, we can approximate |p — k| — > max(|p|, |k|). Thus we write 

j( l L-^,g)=J(l,g)9(k-p) + j(-,g)9( P -k). (3.10) 



max(fc,p) 

Assuming A(p) = A(|p|) and integrating over the angle in Eq. (|3.9[) . we get 



. , x . a f dkkA(k) , 
A(p = A + -j / A2; ^ P, *), 3.11) 

tH y V /A: 2 + A 2 (fc) 



where the kernel 



P \PJ \P 

and K (x) is the complete elliptic integral of the first kind, 9{x) is the Heaviside step function. For zero bare gap, 
Aq = 0, Eq. (|3.1ip admits a nontrivial solution which bifurcates from the trivial one at some critical coupling a = a c . 
To find this critical point we neglect the terms quadratic or higher order in A in Eq. (|3.11|) . It must be emphasized 
that this is not an approximation: it is a precise manner to locate the critical point by applying bifurcation theory 
[37T |. Hence the bifurcation equation amounts to a linearization of Eq. (|3.1ip with respect to the gap function, the 
result reads 

oc 

A(p) = ^ J dkA(k)K(p,k). (3.13) 
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Note that the ultraviolet cutoff, A, has been taken to infinity, which is appropriate at the bifurcation point [37 



This equation is scale invariant and is solved by A(p) 
transcendental equation 



p 7 on the condition that the exponent 7 satisfies the 



4ff 
n 3 Nf 



dx : 



9(1 - x)K(x)J(x, g) + J(l, g) 9 -^—^-K (- 

x \x 



lT 3 Nf 



dx [x -7 J(x, g) + J(l, g)x y - 1 } K(x), < 7 < 1, 



(3.14) 



where J(x, g) is given by Eq. 



and 



2 arccos g 



S<1, 



(3.15) 



Eq. (|3.14[) defines roots 7 for any value of the coupling g. A bifurcation occurs when two roots in (0, 1) become equal. 
Numerically we find that this happens for 7=1/2 and the critical coupling (Nf = 2), 



g c = 1.445, 



(3.16) 



which corresponds to a c = 0.92. For values g > g c the roots become complex indicating that oscillatory behavior of 
the gap function takes over from non-oscillatory one. Equation (|3.14[) determines the critical line in the plane (a, Nf) 
which is presented in Fig[T] This line should be compared with the critical line 



4A n 



(3.17) 



2 - irN f \ c 

obtained in Ref. [23(| using the static polarization function (A c = 1/4 in Ref.[23j] for the kernel approximation 
used below, and A c = 0.23 for more refined bifurcation analysis in Ref. fl7|). The most crucial difference between two 
critical lines is that there is a critical number of flavors, N cr a = 2/ttX c , for the critical line (|3 . 1 7[) for which a = 00 
while a never tends to infinity at finite Nf for the critical line (|3.14[) presented in FigJTJ 

Recently, in Ref. [36| an approximation for the frequency dependent one-loop polarization (|2.5[) was used which 
reduces it to (|2 . 6[) with additional \[2 in the denominator, in this case the critical value a c = 1.13. The more refined 
analysis using bifurcation theory gives a c = 0.93 very close to the value we found above. We remind also that 
renormalization group calculations in two loops yield a c = 0.833 (38[. 



FIG. 1: The critical coupling as a function of Nf. 

A dynamical gap is generated only if a > a c . Since for suspended clean graphene the "fine structure" constant 
a w 2.19 is supercritical, the dynamical gap will be generated making graphene an insulator. Note that for graphene 
on a SiC>2 substrate the dielectric constant k f=a 2.8 and a « 0.78, i.e., the system is in the subcritical regime. The 
value of a c is rather large that implies that a weak coupling approach might be quantitatively inadequate for the 
problem of the gap generation in suspended clean graphene. Therefore, it is instructive to compare our analytical 
results with lattice Monte Carlo studies (25[, where a c = 1.08 ± 0.05 for Nf = 2 that is rather close to our analytical 
findings. 
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IV. NONLINEAR EQUATION AND CRITICAL EXPONENTS 

The above analysis is adequate precisely at the critical coupling, i.e., at the bifurcation point of the original nonlinear 
equation. To study momentum dependence of solutions of Eq. (|3.1ip beyond the critical point we now turn back to 
the case of static polarization when J = n/(l + g) and Eq. (|3.1ip is written in the form 



A 

A(p) = A + ^ J 



dkkA{k) 
v/fc 2 + A 2 (fc) 



K(p,k), A = 



2(l + 7r7V/a/4)' 



(4.1) 



with the kernel (compare with Eq. (|3.12D ) 

AC(p, k) 



1 R f2Vpk 



p + k 



6( P -k) K f 



0(b-p) K (P 
k \k 



(4.2) 



The gap equation is essentially different from a gap equation in BCS theory where the gap is momentum independent. 
In FigI2 we present the results of our numerical solution to Eq. (|4.1[) for A = 0, Nf = 2 and several values of A. The 



A(p) 0.0015 - 




FIG. 2: Momentum dependence of the solution to the gap equation (|4.1|) for Ao = 0, Nf = 2: the bold (black) line for A = 0.3, 
the dashed (red) line for A = 0.27, and the dash-dotted (green) line for A = 0.25. 

gap is weakly dependent on a momentum up to values p ~ A(0) after which the behavior becomes more steep. To 
estimate the gap A(0) we need to know the bandwidth parameter A which can be obtained by equating the wave vector 
integral over the Brillouin zone to the integral over two Dirac points with a cutoff at k c . We get k c = (7r/v / 3) 1 ^ 2 (2/a) 

where a is the lattice constant, therefore, restoring H and vp = s/3ta/(2h), we find A = hvpk c = \J 7V\/3t ~ 2.33t. 
For the hopping parameter t = 3cV, we obtain A w 7eV. The maximal possible gap A(0) is reached for a — > oo 
that corresponds to the value A = 1/tt sa 0.32 (Nf = 2). For the values of A's used in Fig. [2] we find the estimates: 
A(0) = 200K,40K,25K for A = 0.3,0.27,0.25, respectively. 

To get insight into analytical solutions of Eq. (|4.ip we approximate the elliptic integral functions in (|4.2p by their 
asymptotical values at p -C k and p 3> k, we obtain the kernel 



7T (6(p-k) 6(k-p) 



P 



(4.3) 



This allows us to reduce the nonlinear integral equation (|4. 1 [) to the second order nonlinear differential equation 



with the infrared (IR) and ultraviolet (UV) boundary conditions: 



p 2 A'(p) 
( P A(p))' 



p=0 



= 0, 



= A . 



(4.4) 

(4.5) 
(4.6) 
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Eq. (|4.4[) is scale invariant, i.e., if A(p) is a solution, then lA(p/l) is also a solution. The scale invariance is broken by 
the UV boundary condition only. 

The chiral condensate (0|4 n I / |0) is the order parameter for the semimetal-insulator transition in graphene, it breaks 
spontaneously the initial U(2Nf) symmetry down the U(Nf) <g> U(Nf) but keeps parity and time-reversal invariances. 
It is expressed through the full fermion propagator as follows 



1**10) = -to<0|T*(a;)*(0)|0) = -itr J g j j^G(uj,p) 

(4.7) 



7T J y/p 2 + A 2 (p) 7rA 



p=A 



where for the last equality we used Eq. (|4.4p . Hence the condensate is nontrivial if a nontrivial solution of the gap 
equation exists. 

One can easily find the solutions of Eq. (|4.4[) in two asymptotic regions. For p <C A(p), 

A(p)=Ci + ^. (4.8) 
P 

The IR boundary condition (|4. 5|) implies C2 = 0, therefore, A(p) ~ C\ for p <C A(p). For p ^> A(p), 

A(p)^C^-T++C7 4 p-T- ) 7 ± = ^±Va^aT. (4.9) 

Clearly, in order to find a solution of Eq. (|3.11|) one needs to show that there exists a solution of the nonlinear differential 
equation (|4.4|) which connects the asymptotic A(p) ~ const in the infrared region, p — > 0, with the asymptotic (|4.9j) 
at large momenta. For this, let us define 

A(p) = e*u(i + f ), t = lnp, (4.10) 
then the function it(i) satisfies the differential equation 

u" + 3^ + 2m + X - = (4.11) 

The IR boundary condition implies 

e 2t (u' + u) =0. (4.12) 

i— — 00 

We require that e t u(t) — » 1 as t — > —00 since all other solutions for A(p) are obtained by varying the constant <o- For 
this normalization, the infrared scale for the general solution is given by A(0) = e~ ta . 

The dependence of the integral equation (|3.11[) on the bare gap Ao now becomes an ultraviolet boundary condition 
for the differential equation, it is 

u'(t A + t ) + 2u(t A + t ) =A /A. (4.13) 

This condition determines the value of the parameter to = — lnA(0) as a function of the coupling constant, A, the 
bare gap, Aq, and the cutoff, A. Eq. (|4.11j ) can be rewritten in the form 



or, equivalently, 



+ 3u' = — 7-V(u), V(u) =w 2 + A v / l + u 2 . (4.14) 
du 



hu'Y + V(u)) =-3( U ') 2 - (4.15) 



Eq. ()4.14p is the equation for a particle of unit mass moving in a potential V with friction proportional to velocity. 
The "energy" \{u') 2 + V(u) reaches its absolute minimum at u = 0, hence the particle moves toward u = damped 
by the friction. The asymptotical behavior near u — is described by the linearized equation 

u" + 3u' + (2 + X)u = 0, (4.16) 
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and depends on whether the coupling A > A c = 1/4, or A < A c 

B 



u(t) -> -=== e 3t/2 sinh y/ \ c - A (t + 5) , i — > oo, weak coupling (A < A c ), (4-17) 
v A c — A L J 

u(t) — > e~ 3t/2 sin V A - A c (t + S) , i — * oo, strong coupling A > A c , (4-18) 
V A — A c L J 



where the constants A, B, and <$ are functions of the coupling constant A. We explicitly singled out the factor l/\/A c — A 
in front of Eqs. (|4.17p , (|4.18p since the function u(t) must be nontrivial at A = A c . Obviously, A(X = A c ) = B{\ = A c ). 

The asymptotics (|4. 17[) and (|4.18[) imply that at weak coupling the particle situated initially at u(— oo) reaches 
u = for infinite time, meanwhile, at strong coupling it will get to u = in a finite time and then oscillate there with 
damped amplitude. 

It is easy to see that at weak coupling there are no nontrivial solutions satisfying the UV boundary condition for 
Ao = 0. For strong coupling, the UV boundary condition (|4.13[) with Ao = admits an infinite number of solutions 
for the gap scale A(0), corresponding to different solutions of the equation 

A\/\ 

u' (t A + h) + 2u (t A + to) « , e -3(*A+to)/a sin [e + < f) = 0, (4.19) 

V A — A c 

where 

9 = s/\ - A c (f A + t + 5) = VA - A c In (^y) > 4 = arctan 
Hence, the solution is given by 9 = nn — <fr, or 



(4.20) 



A(0)^A e "cxp , n=l,2,.... (4.21) 

The solution without nodes, n = 1, corresponds to the ground state since it generates the largest fermion gap and 
has the lowest energy. The critical coupling A c = 1/4 is a bifurcation point of the integral equation (|3.1ip with the 
static vacuum polarization [39l ]. The expression (|4.2ip for the gap implies that this bifurcation point corresponds to 
a continuous phase transition of infinite order. As was shown in Ref. [T4I, the critical coupling A c is closely related to 
the phenomenon "fall into the center" in quantum mechanics problem. A similar situation takes place in the strong 
coupling QED4 0] where in the ladder approximation (and more generally in quenched approximation when fermions 
loops are neglected (40j |) the phase transition is also of infinite order. The dimensionless correlation length, 

(-m-"hrx)- (4 - 22) 

exponentially grows when A — > A c . Such a behavior is inherent for the Berezinskii-Kosterlitz-Thouless phase transition 
(or the conformal phase transition (CPT) (26|) and, obviously, is related to the scale invariance of the problem under 
consideration. Note, however, that taking into account the finite size of graphene samples should turn this phase 
transition into a second order one (as it is shown for QED3 in Ref . [28| ) . 

Eq. (|4.22p means also that in nonperturbative phase there is a nontrivial running of the coupling a (or the Fermi 
velocity vf) though we neglected its perturbativc running. Defining the f3 function in a standard way, we find 

P(a)^A^ = ~(l + nN f a/4) 2 (\-\ c f/ 2 , A > A c , (4.23) 

where A is defined in Eq. (|4.ip . The f3 function depends nonanalytically on the coupling a and can not be obtained in 
perturbation theory. We expect that if a pcrturbative running of a is included the critical point A c becomes a second 
order phase transition point on which the /3 function (|4.23p is continuous when approached from both perturbativc 
and nonperturbative phases. 

The order parameter (^^>) in terms of the function u(t) is given by 

= -± e 2t * («'(i A + t ) + u(U + h)) , (4.24) 

7rA 



and equals 



= »1 A A^A^Wsin^) = -^AVW/^O), (4.25) 

tt v /A(A- A c ) ttA 3 / 2 



9 



where the relation 9 = tt — (f> was used. 

For nonzero bare gap, A 7^ 0, we obtain the following equation for the scale A(0): 

aV\ a 3 / 2 (o) . ,„ . „ 

A = , , sin 1 



and the order parameter, 



— N f 

7TA 



An- 



A A 3 / 2 (0) 



sin ( 



(4.26) 



(4.27) 



Let us write 8 
as 



V / A^A C " VA 

tt — e where e tends to zero when Aq — > and A — > A c . Then the above equations are rewritten 



An = 



AV\ A 3 / 2 (0) 



VX^Xc \/A 
_ N f A \2X- 1 



-A 



2A 



sine, 



.4 



A 3 / 2 (0) cos< 



VA Va 



(4.28) 
(4.29) 



In such a form the equations are convenient for finding critical exponents near the phase transition point A c . Critical 
exponents describe the approach to criticality of such quantities as the correlation length, the order parameter, the 
susceptibility, etc., they are defined in a standard way [Tel l4l|: 



A 



A(0) 



(A-Ao)- 



(#40 



A 



A=A C 



l/S 
> 



A 



* 0. 



(A — A c ) , X 



dA Q 



(A - A c ) _ 



A — > A c , 



(4.30) 
(4.31) 



If the theory of second order phase transition is applicable, then the exponents are assumed to obey the following 
hyperscaling relations in spaces of arbitrary dimension D: 



2(3 + 7 = Dv, 2/35 - 7 = Dv, 



8-1 2-Tj 



5+1 D 

Here the exponent 77 describes the behavior of the correlation function 



(3 = v 



D-2 + ri 



(##(r)##(0)) 



-D+2-r] 



Using Eq.g^U), we find 



and due to Eq. (|4~26"l) . 



A=A C 



4iV/A 



A=A C 



An 



2AA 3 / 2 (0) 



A(0) 



An 



2/3 



A = Ac 



>(A/A ) 

the critical exponent 5 = 1, and from hyperscaling relations we obtain 

»7 = 2, 7 = 0, /?=y. 

The infinite order phase transition with the correlation length (|4.22j) formally corresponds to the limit 



00. 



(4.32) 



(4.33) 



(4.34) 



(4.35) 



(4.36) 



(4.37) 



Certainly, the infinite order phase transition is quite different from that one studied in lattice simulations |25j | where 
a second order phase transition was found with the critical exponents 5 ~ 2.3, (3 ~ 0.8,7 ~ ^-V^S = 2). One of the 
reasons for such a difference might be a finite size of a lattice that changes the kind of phase transition. Effectively, 
the finite size of a lattice can be taken into account by introducing an infrared cutoff (kg ~ tt/L, L is a linear size 
of the system) in the integral equation (|3. 1 1|) [H, ■ Another reason could be that one should take into account 
residual lattice interactions, i.e., the present analysis has to be further refined by incorporating effective four-fcrmion 
terms (see the next section). 
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V. PHASE DIAGRAM IN THE MODEL WITH ADDITIONAL FOUR-FERMION INTERACTION 



As discussed in Introduction, when comparing the results of lattice simulations [2a ] with analytical calculations 
one should have in mind that the continuum theory described by Lagrangian (|2.ip putted on a lattice contains 
unavoidably additional interaction terms, in particular, local four-fermion interaction terms. This means that it 
would be appropriate to add to the continuum theory some local four-fermion interaction terms in addition to the 
long range Coulomb interaction in Eq. (|2.ip . The amount of induced couplings depends of course on the particular 
lattice regularization employed. Furthermore, according to (3ll . I32I f33j ] , the effective continuum model for quasiparticles 
in graphene in addition to the Coulomb interaction should contain contact four-fermion interaction terms which arise 
from the original lattice tight-binding model. Usually these residual four-fermion terms are irrelevant operators from 
the point of view of the renormalization group, however, we will show that they can play a significant role in the 
critical behavior. In order to study how these four-fermion terms influence the gap generation, we will consider in this 
section a continuum model with the Coulomb interaction and the Gross- Neveu four-fermion interaction of the type 



C 4 



f (*a*«) 2 , 



(5.1) 



where the four-fermion coupling constant G is of the order of the lattice constant and the "flavor" index a = 
1,2, Nf, Nf = 2 for physical spin-1/2 electrons. The interaction term (|5.1[) breaks the initial U{2Nf) sym- 
metry of the action (|2.ip down to the U(Nf) ® U(Nf) (g> Z^ symmetry. While the gap term A'J"!' is invariant under 
the U(Nf) <g> U(Nf), it is not under the discrete chiral Z2 - symmetry: VP — > 75 'J, \& — > —^75. In the absence of the 
bare gap term, Ao x I n I / , the Zi symmetry forbids the fermion gap generation in perturbation theory. The appearance 
of the energy gap is due to the spontaneous breaking of the above discrete chiral symmetry that leads to a neutral 
condensate (\I>\I/) of fermion-antifcrmion pairs (excitonic condensate). 

The gap equation (|3.1ip is modified in the presence of the interaction (|5.ip in the following way: 



A(p) = A - G(l 



1 



AN 



■)<**) 



a f dkkA(k) 



^k 2 + A 2 (k) 



JC(p, k), 



(5.2) 



where the condensate contributes like a bare fermion gap and can be computed from the fermion self-energy; 

the factor 1 — 1/4JV/ in the second term on the right hand side takes into account both Hartree and Fock (—l/4Nf) 
contributions. For the sake of simplicity we consider only the Hartree term, if necessary the Fock contribution can 
be easily restored in final formulas [1^. In the approximation to the kernel used above (Eq. (|4.3p ) the condensate is 
given by the expression (|4.7p . The condensate does not change the differential equation (14. 4p . however, it modifies 
the ultraviolet boundary condition (|4.6[) : 



1 



gNf 



pA'(p) + A(p) 



= Ar 



(5.3) 



where we introduced the notation g = GA/tt and A is defined in Eq. (|4.4p . Using the definition of the gap function 
in terms of the u(t) function (|4.10p and the asymptotic behavior of the last one, Eqs. (|4.17p . (|4.18p . the equation 
can be written for A 3> A(0) in the following form: 



B 



D 



A 



A 3 / 2 (0) 



A 

A 3 / 2 (0) 

n/A 
A 3 / 2 (0) 



A 



1 



gNf 
X 

qNl 
X 

gN f 



cosh ijj In 



+ 



1 



Ae 5 \ 

gNf IX Ae s 1 
In 



2w V A (°) / 



A , w = vAc 



~ A(0) 

cos u> ln 



= An, A = A, 



A, A < A c> (5.4) 
(5.5) 



Ae s \ 



A(0)J 



1-gNf X . 

1 sin oj ln 

2Q 



Ae s \ 
A(0)J 



A , 



\/X-X c , X > A c , (5.6) 



and we remind that in the utilized approximation A c = 1/4. These equations imply the following solutions for the 
dynamical gap in the case An = 0: 



A(0) = Ae s 



gNf(l-2u)-\(l+2w) 



A(0) = Ae" exp 



A(0) = Ae-'exp 



SfiV>(l + 2uj) - A(l - 2ui) 
~gN f + 1/4 



A < A c , 



-2 



gN f - 1/4 
7m 1 
Co oj 



A = A c = i , 



n ~gN f + X 

arctan 2uj 



gNf - x 



gN f > i 

, A > Ac, 



1,2, 



(5.7) 
(5.8) 
(5.9) 
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Note that the solutions (|5.8p . (|5.9[) contain essential singularity, the first solution at g = 1/ANf and the second one 
at A = \ c . 

Setting A(0) = 0, we find the critical line separating the spontaneously broken and unbroken phases of the chiral 
symmetry: 



9 = ^(i + Vi-VAc) 2 , for A<A C = ±, 



9 < TNji ^ — A c 



(5.10) 



The phase diagram in the plane of two coupling constants is displayed in Fig[3] [42j]. Above the critical line the gap 

0.5 4 



**) ^ 




FIG. 3: Phase diagram for N f = 2. 

equation for the fermion self-energy A(p) has a nontrivial solution. Thus, the chiral symmetry is dynamically broken 
that implies the existence of a nonzero vacuum condensate ('I"!'). For g = 0, the condition for a gap generation 
becomes A > A c and the corresponding critical coupling coincides with Eq. p.l7[) . On the other hand, in the other 
limiting case a = 0, g c = 1/Nf coincides with the critical coupling in the Gross-Neveu model. In the part of the 
phase diagram above the critical line A < A c the short-range four-fermion interactions play important role for the 
condensate formation, meanwhile, in the region A > A c Coulomb forces are mainly responsible for the condensate 
formation. 

We consider now the phase transition along the upper part of the critical line and compute the critical exponents. 
Since we consider non-running coupling a (in absence of running of the Fermi velocity vp) the renormalization group 
flow can be determined from Eq. (|5.7p which near a critical point takes the form 



A(0) 
A 



9-9i 
9-h 



9i 



AN f 



92 



(1 - 2uf 



9 > 9i > 92- 



It implies an explicit form of the (3 function for the coupling g: 

= -Nf(g-gi)(jj-gn), g > gi- 



q,A(0) 



(5.11) 



(5.12) 



Eq. (|5.12[) has indeed a nontrivial fixed line at g = g\ . We stress that the (3 function (|5.12[) is obtained in nonperturba- 
tive phase where a quasiparticle gap is spontaneously generated. In perturbative phase the j3 function was calculated 



in [43| (see, Eq.(7) there), in the leading order in 1/Nf and small coupling a both (3 functions behave as (3 ~ 
near the fixed point go ~ (1 — a)/Nf. As is seen from Eq. (|5.11[) . the phase transition is of the second order, 
the deviation from the critical line as r = 



gi and because A(0) ~ r ' (r — > 0) we find the exponent 
1 



The condensate is given by 



NfB 
ttA 



A 3/2 (0)A l/2 



cosh lo In 



Ae 5 



sinh 

2ui 



i In 



Ae s 
A(0) 



-(.9 -go) 

Denoting 



(5.13) 



(5.14) 
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On the critical line, Eg.([5"3]) implies A(0) ~ Ao /(3/2+w) . Substituting it into the expression for the fermion condensate 
gives the critical scaling relation 

3/2 + w 



thus the critical exponent 



3/2- w' 



(5.16) 



It is equal to <5 = 2 in the case of the Gross-Neveu model when a = and 5 — ► 1 for a — » a c . It is easy to find also 
the /3 exponent: 

Finally, it follows from Eqs. (f5T4|) . (|5~T4|) that 

Sa <**) r^O, (5.18) 

A o =0 

hence the exponent 7=1. The found critical exponents satisfy the hyperscaling relations (|4.32|) . The additional 
critical exponent rj may be calculated independently, or using hyperscaling relations, r] = 2 — 2u>. By definition, the 
anomalous dimension j rn of the composite operator is given by dim(\E r \E r ) = D — 1 — 7 m , then the correlator (|4.33[) 
implies the relation r\ = D — 2^ m . In our case, D = 3, we obtain 7 m = 1/2 + w. The dynamical dimension of the 
four-fermion interaction term equals dim(>I"I') 2 = 2dim (\f' 1 I r ) = 4 — 27 m . Because 1/2 < 7 m < 1 along the critical 
line, dim(\f r \I/) 2 < 3 and the four-fermion operator (>I"J') 2 acts as a renormalizable one. In the renormalization group 
terminology, the (\E"J') 2 becomes a relevant operator in the scaling region while it is irrelevant away from the critical 



line, in accordance with standard renormalization group approach 32j, as its effects are suppressed by powers of 
cutoff. On the other hand, the anomalous dimension "f m governs the behavior of the amputated Bethe-Salpeter wave 
function (form factor) of bound states, ^( am P) (q) ~ ((7/A(0)) 7m_1 , in the range of momenta A(0) « 9 « A. The 
"critical" value -f m = 1/2 separates loose (-f m < 1/2) and tight (j m > 1/2) bound states. The wave functions with 
large j m (-f m > 1/2) slowly decrease with momentum, they describe tight bound states which are relevant for critical 
scaling laws of a theory [44j |. Since such bound states resemble point-like particles, the scaling properties of a theory 
can be described by an effective Lagrangian with elementary scalar fields (for recent such an approach, see, Ref . fill] ) . 
The computer simulations of lattice graphene model may reveal in principle the existence of such tight bound states. 

We see that the additional Gross-Neveu four-fermion interaction plays an important role: First, it changes the 
order of a phase transition from the infinite to the second order one. Second, the critical coupling becomes lower 
than in the model with the pure Coulomb interaction. Third, the critical exponents stay closer to those obtained in 
lattice simulations 25]. Further, the critical indices depend on the coupling a along the critical line < a < a c and 
satisfy the hyperscaling relations. The phase diagram (|5. 10[) resembles closely those obtained in the strong coupling 
QED4 [4(| and QED3 (47[- Since the phase transition is of second order along the < a < a c part of the critical 
curve, Eq. (|5.10p . resonances should exist on the symmetric side of the curve, whose masses tend to zero as the critical 
curve is approached [H, Eii [.The part of the critical curve with g < 1 /4Nf is rather special and is related to the 
conformal phase transition 26j . It is characterized by a gap function having an essential singularity at the transition 
point, and by abrupt change of the spectrum of light excitations as the critical point is crossed: light bound states 
near the critical line are absent in the symmetric phase, however, they are present in the phase with broken symmetry 
(for a discussion in detail of the CPT in QED3, see Ref. [48l [50j] ) . The corresponding effective potential for the order 
parameter (^W), unlike the familiar Ginzburg-Landau potential, was shown to have a branched fractal structure in 
the region a > a c , where the Coulomb interaction is mainly responsible for the bound states formation [24| . 

We hope that at least some features of the picture outlined above will be confirmed in experiments with suspended 
clean graphene. 



VI. CONCLUSIONS 



In this paper we studied the gap generation in suspended clean graphene at neutral point. Solving the Schwinger- 
Dyson equation with the frequency dependent polarization function we found analytically that the critical coupling 
constant for onset of a gap equals a c = 0.92 which is close to the value obtained in Monte Carlo simulations. We 
showed that the critical coupling a c corresponds to the infinite order phase transition in the case of purely Coulomb 



13 



interaction with peculiar critical exponents while Monte Carlo simulations point to the second order phase transition 
with different critical exponents. 

Adding the Gross-Ncveu four-fermion interaction that is present in the continuum limit of the lattice model we 
found the critical line Eq. (|5.10p in the plane of Coulomb and four-fermion coupling constants separating zero-gap and 
gapped phases. We showed that the order of a phase transition changes from the infinite to the second order one 
along the part < a < a c of the critical line and the critical coupling becomes lower than in the model with pure 
Coulomb interaction. The critical exponents v, <5, (3 along the line of second order phase transition are given by the 
expressions (|5.13[) . (|5.16[) . (|5.17p . respectively, and the exponent 7=1. These exponents satisfy hyperscaling relations 
and characterize the transition between phases with distinct symmetry properties and become, in general, functions 
of the Coulomb coupling a, or the four-fermion coupling g. They are close to the critical exponents obtained in lattice 
simulations p5j . 

The other part of the critical curve with g < 1/ANf is rather special and is related to the conformal phase 
transition characterized by an essential singularity at the transition point and by abrupt change of the spectrum 
of light excitations as the critical point is crossed. However, the shape of the last par t of phase transition curve 
might be strongly influenced by the finite size effects which appear to be nontrivial [281 ] - Also, the running of the 
coupling a, due to the running of the Fermi velocity vp, may change the shape of the vertical part of the critical 
line. These effects most likely change the kind of phase transition to the second order one. This would indicate that 
the semimetal-insulator transition in graphene is likely to be of second order. We expect that the form of the critical 
curve in graphene can be checked in further lattice simulations. 

Also, our results maybe important for the proper interpretation of lattice simulations of low-energy field-theory 
model for quasiparticles in graphene interacting through the long-range Coulomb potential [25| because local four- 
fermion terms are expected to be generated by the lattice rcgularization procedure. We showed that in spite being 
small (G ~ lattice constant) the induced local interactions can play a significant role in the critical behavior observed 
in lattice simulations. A related aspect of the near-critical behavior is the appearance of composite electron-hole 
degrees of freedom whose form factors slowly decrease with momentum (tight bound states), and the momentum 
behavior is governed by large anomalous dimension. Their dynamics can be studied similarly to that in strongly 
coupled QED [5l[ but this remains a problem for future investigations. 
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